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Abstract 


A model that yields the spatial correlation structure of 
atmospheric mass field forecast errors has been developed. The 
model is governed by the potential vorticity equation forced by 
random noise: 

(v2-Co sin^O) ^ (X,0;u) = F(X,0;o)) (1) 

where v2 is the Laplacian operator on the unit sphere , A and 9 are 
longitude and latitude, ^ is the geopotential error field at 500mb 
and P is white noise corresponding to a random realization «. 

The spatial covariance function r is defined by 

r ’CXi , 0; A2»e2) “ E (Xi,0i;a)) ^ (A2,02;w)}, (2) 
where E is the expected value. 

Three methods of solution have been tested. In the first 
method, Eq. (1) was solved by expansion in spherical harmonics 
and the correlation function was computed analytically using the 
expansion coefficients. In the second method, the finite-dif- 
ference equivalent of Eq. (1) was solved using a Fast Poisson 
Solver. The correlation function was computed using stratified 
sampling of the individual realizations of F(w) and hence of 
(j) (o)). In the third method, a higher-order equation for r was 
derived from Eq. (1) and solved directly in finite differences 
by two successive applications of the Fast Poisson Solver. The 
three methods were compared for accuracy and efficiency, and the 
third method was chosen as clearly superior. 

The results agree well with the latitude dependence of ob- 
served atmospheric correlation data. The value of the parameter 
Co which gives the best fit to the data is close to the value 
expected from dynamical considerations. These results provide 

the basis for an optimal choice of coefficients for statistical 
analysis of atmospheric data. 
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1- Introduction 

The statistical structure of large-scale atmospheric 

fields is of both theoretical and practical interest to meteor- 
¥ 

ologists. Theoretically, it is of interest to know what this 
structure is and how it becomes established (Gandin, 1963). In 
particular, the connection between the atmosphere's dynamics and 
its statistics is an attractive area of study. 

Practically, numerical weather prediction (NWP) requires 
the accurate, detailed description of atmospheric fields as a 
starting point for their forecasting. The data available for 
such description are nonuniformly distributed in space and 
contaminated by various errors {Bengtsson, 1975). It is necessary, 
therefore, to use some form of interpolation to derive field values 
at the points of a uniform grid. It is desirable, furthermore, 
that these values be as free of errors as possible. 

Interpolation coefficients can be chosen which will minimize, 
under certain assumptions, the expected value of the interpolation 
error, given the statistical properties of the errors in the data 
(Rutherford, 1972). This statistical approach to meteorological 
interpolation has become increasingly attractive recently, due 
to the large number of different data sources with varying error 
characteristics made available by the Global Atmospheric Research 
Program (GARP) (Fleming et al . , 1979). It is often referred to 
as "optimal interpolation" (01) and has been implemented opera- 
tionally by the O.S. National Meteorological Center (NMC; 

McPherson et al., 1979), among others. 
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The optimal choice of interpolation coefficients ir 01 
clearly depends on knowing the statistical properties of the 
fields one wishes to interpolate. Hence the practical impor- 
tance of an accurate model for large-scale atmospheric statistics. 
The purpose of this report is to contribute bo the formulation 
and validation of a dynamically based model for atmospheric 
statistics. 

In a study on the use of satellite-derived temperatures for 
NWP, Ghil et al. (1979; to be referred to as GHA) were led to 
consider the difference between the observed atmospheric tem- 
peratures, T°, and model-forecast temperatures, T^, The model 
used in that study was the nine-level 4" lat. x 5“ lone,, primi- 
tive equation, spatially second-order model of the laboratory 
for Atmospheric Sciences of NASA’s Goddard Space Plight Center 
(GLAS); temperature data were obtained from the Data System Test 
DST-6 held during January-March 1976. 

The spatial correlations of the difference field T*^ - T^ 
were computed. It turned out that, for the same spherical distances 
s between points, correlations were typically higher in the 
tropics than in high latitudes. In other words, the correlation 
r(g,n) of temperatures T° - T^ at a point 5 on the Earth with 
those at a point n a distance s away, s = disb (C,n), falls 
off more rapidly with s the higher the latitude of the point r, 

(Pigs. 2a-d, GHA). No large or systematic dependence on height 
was observed when stratifying the correlations by pressure level 
rather than by latitude. 


„ . , t:. . 
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It was suggested ( ibid. ) that this striking latitude depend- 
ence of the studied field's second-order moments reflects the 
dependence on latitude of the Rossby radius of deformation, L. 

''The latter is a characteristic length scale for a number of dynamic 
phenomena which determine the spatial structure of atmospheric 
fields. Wo decided to pursue this heuristic suggestion further, 
and formulated the stochastic-dynamic model investigated in this 
report. 

Section 2 presents the model, and the governing equation, 

Eq. (2.7). This equation is solved by a series expansion in 

Section 3. For given, fixed right-hand side, Eq. (2.7) can be 

solved numerically by the use of a generalized Past Poisson 

solver, as shown in Section 4. The full stochastic form of 

(2.7) is solved by Monte-Carlo simulation in Section 5. An 

equation for the covariance funtion solution 

* 

to (2.7) is derived in Section 6. It is seen to depend on a 
scale parameter, c«, 

^ r = r(?3_, ^ 2 ' • 

Comparing model correlations with observed mass field correlations, 
we obtained the best value of Cq. Numerical results are presented 
in Section 7. Concluding remarks follow in Section 8. 
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2. DynamicaX model of the forecast error field 

We will assume that for periods of a few days, the dynamics 
of the atmosphere are approKimately governed by the equation of 
conservation of potential vorticity 

d ( a“2 V 2 ^ ^ 4 , f)«Q, (2.1) 

dt 

Here t}i « 4> /£q is the quasi-geostrophic stream function, ^ being 
the height h of the 500mb surface multiplied by the gravity 
g, «l> *gh; the Coriolis parameter is f«2 n sin e , with n the 
angular velocity of rotation of the earth and $ is latitude, 
while fo is a constant value of f corresponding to a mid-la titud 
e 0Q. The radius of the earth is a, 7 ^ is the laplacian 
Operator on the unit sphere, with X longitude, L is the Rossby 
radius of deformation, L 2 =gD/f 2 , with D a characteristic depth. 

The forcing term Q represents diabetic heating, dissipation and 
lower-boundary effects. 

Equation (2.1), with Q=0, is strictly valid for a quasi-geos- 
trophic, frictionless, shallow-water model without topography, 
with a mean depth D. It is also valid for each of the vertical 
modes in a linearized quasigeostrophic model in which the vertical 
dependence has been separated out (Phillips, 1973). In this 
case D is the equivalent depth corresponding to either the external 
mode or to one of the internal modes. Our assumption is that 
for periods of a few days equation ( 2 . 1 ) is a reasonable model of 
large-scale atmospheric flow. 

The dynamics of a numerical weather prediction model are also 
governed by an approximation of equation ( 2 . 1 ): 

rs# rw 

d (a "2 7 2 -l -2 ^ +f) = ^ 2 . 2 ) 




rftf I 'fi i n* 




§. 

I 

The tilde represents the numerical truncation effect in the operators 
d/dt and V 2 on the one hand, and the errors in the parameterization 

IV 

of the physical forcing, Q, on the other. 

It follows that forecast errors 6\|) ■ will also be 

(governed to a good approximation by a conservation equation of 
potential vorticlty, which does not contain the planetary vorticity 
term f : 

d ( a~2 -L"2 6ij;) «sorrors; (2.3) 

n't 

we let the errors in the right-hand side of (2.3) represent all 
the approximations, physical and numerical, made in equation 
(2.1) and, a fortiori, in equation (2.2). At the initial time, 
t«0, ^ is obtained from observations of the atmospheric state 
which are also made with certain errors: 

=srrors at t“0. (2.4) 

If the errors in both Eq. (2.3) and the initial conditions (2.4) 
were zero, then one would obtain that the potential vorticity of 
the error will remain identically zero 

(a-2 v2 -l-2) =0 

at all times t > 0. 

In the presence of purely random errors , we can combine 
Eqs. (2.3) and (2.4) to yield a time-independent equation governing 
forecast errors. This equation is a stochastically forced steady- 
state potential vorticity equation on an f-plane: 

(a-2 7^ -L”2) 0 , X; w ). (2.5) 

We take F-t to be random white noise, corresponding to different 
realizations of atmospheric processes labeled by uD,at time 
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t>Q, We expect Bq. (2,5) to be a good representation of the 
structure of the forecast error field whenever systematic errors 
in Eq. (2,1) and (2,2) are small. It will cease to be valid at 
time scales longer than a day or two| because the errors in the 

JV 

approximation of the nonlinear terms, d/dt-d/dt, as well as in the 
parameterization of physical processes, Q-Q, become sizeable and 
nonrandom. We may also expect that Eq, (2,5) will be less accu- 
rate in the Northern Hemisphere than in the Southern Hemisphere, 
because topographic forcing and land-sea contrast are more impor- 
tant in the former than in the latter. 

The statistical properties of Ft, in particular its variance 
might change with time. Since t in (2.5) is only a para- 
meter. we shall consider a fixed in the sequel. The value 
of a affects only the amplitude, and not the structure of the 
solution. 

The Rossby radius of deformation L depends on the equivalent 
depth D and on the sine of latitude; 

L-2s (4 n2/gD) sin^ 9. (2.6) 

For simplicity, we assume that one vertical model dominates the 
error field, and shall determine the value of the equivalent 
depth that best fits the data. 

Summarizing, we will study the equation for the geopotential 
error field 4>, 

(V2 -Cq sin2 0 ) ({. =F(X,0jw), (2.7) 

where Eq. (2,5) was multiplied by the constant f^a^, and P is a 
spatially multi-dimensional noise process; 

e{p(€; to)l-0-, (2.8a) 
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E{F(Cj^;w) P (2.8b) 

Here E io tho expectation operator or ensemble average over the 
individual realization », is the position vector, and 

o2 a prescribed variance. 

We are mainly interested in the covariance function. 

P(5i, 52) - E (KClSw) ^ (52SW)} (2,9) 

of the solution ^ (€,‘w). The reason for our interest in r is 
that intei’polation formulae for assimilation of atmospheric data 
require r es the basic statistical information. 

In the following four sections, our methods for the solution 
of (2.7) and for the (.computation of (2.9) are described and com^ 
pared . 
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3'. Sggiss lolut^lon of the model equation 


Consider expansion of , 

MBA) - I z b; y; , 

n»0 m**-n 

in thn Bpherioal harmonies , 


C3sl) 


vJC0/X) « P^(sin0) - 0^^ P™(li) . (3,2a) 


Here y « sin0 and P^(p) are the associated I^egeridre functions 
(Courant and Hilbert, li)53? Hobson, 1955), 



u 


2 . m/2 
2%i 


_ti4*yn 

ap"+™ 




(3.2b) 


normalized so that 



4 7r (2n^-l) (n~m) 1 
(h+m) 1 


'^n,n' 


A 


m,m 


5 I 
n,n' 


(3.3) 


* 

{ ) denotes complex conjugation , 

The representation (3.1) would diagonalize Set. (2.2), provided 

the operator in (2.2)were a pure Helmholtz operator, i.e., only a 

2 

constant, 0 -independent term were added to the Laplacian v . As 
it is, we shall show that (3.1) leads to a five-diagonal representation 
of (2.2). 


orthogonal functiona P^(p) satisfy the three-term recursion 
relation 

'''’n ■ 2^1 *'n-l 5H?r '’n+1 ' £o>^ " > « • 


It follows that, for n > 1 , 


p“ - 

S® pS + 

%-2 ' n-2 ^ 

a® pS 

^n+2 ^H®.2 ' 

(3.4a) 

where 

a® „ » 1 
n-2 

^ n+s\ 
^2n+l/ 

/n+s-1 1 
^2n-l /' 

(3.4b) 


3® - 1 

^n+s \ 

• 

/ n-s \ / n-s+l\ /n+s+1 ^ 

, (3.4c) 


^2n-fy 

^2n-iy ^2n+1^^2n+3 y 

and 

Y® » i 

’n+2 1 

f n-s+1 

2n+l 
i ' 

\ fn-s+2 ^ 

/ ^2n+3 J ' 

(3.4d) 

Note that 

**Q*" 2 
ft ^ » 

q 

s"g-i , 

q 

0 1 

“ Yg ** Yg * 0 for 

q ^ 0. 


substituting the series (3.1), truncated at n»N#into 
Eq. (2,2) we have 


(7^ - c- sin^9) E 


N n 


(V^ - 


(V^ - 


Co e) z z ®n ''n ' 

n«Q m--n 


p„v^) [b® V® + sf y-1 + B® V® + b\ yj 


N n _ -1 

+ y I y"' == F. 

n=2 m=“n 


(3.5) 


i‘rrnn wifrr irrr'illiirrai^-*^ - 






■iiitrir'-iyMMii' 
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On the unit sphere, 


yJJ - - n(nA)yJJ . 


(3,6) 


The fi'.*. , few terms of the five-term recursion (3.4) are given 
by 


pO ^ ^2 ^ I 


I (3y^-l) 


•f ^ ss ~ p2 + ~ p® bs y® -f i. y® 

3 3 *^2 ^ 3 0 3 ^2 ^ 3 0 ' 


yJ = p^ pj “ « f 


|(5 p^-3p) 


X = 2 -0 . 3 n(J « 2 v(J X 3 vO 

+ 5^^ “ I P 3 + ^ ^1 *“ 5 ^3 + “5‘ ' 


P^ ^1|^” e"^^p^ p“^= e"^^pp“^= e"^^ 


i p“^*+ i. p"^ 
5 ^ 5 ^3 




and 

p^ y^; = p^P^ = ^^^p 

Thus (3.5) can be rev;ritten as 
0 


- • 

1 ni 

1 ^il 

*■ m 

3 „X_, 2 „1 

3 2 

0 • 

= 3 ® 

? ^1* 5 '■3 


I Y'“lx i v-1 
5 1 5 *3 ' 


3 1,2 1 


B 


[- °0<i'’'2 + + B-i[- 




] 


[-2 - Co(| Y° + I Y°)] + -2 y^ - CQ(f5 yJ + 1-5 yi)] 

m=-n ‘^n-2]^n-2 [ 'n(n+l) “ Cq 6^ 1 Y^ + 

■^ [ " ^ 0 ^n+2] ^n+2 j 


P . (3.7) 


u. 


Define 


m 

n 


a. 


.m 


n-2 


[-n(„«) - i;; ]aU , 

[■ °0 “n-2 j*n-2 ’ 


(3. 8a) 
(3.8b) 


and 


ra 

'n+2 


- C Y™ T a”' 

I 0 ’n+2 ]^n+2 


(3.8c) 


. * 

Multiplying both sides of Eq.(3.7) by (Yj^) and integrating, we have 

;2ir fl . ^ f2tr /I . * 

dX j ^(yJ) {LHS} dy = dX j ^ (yJ) P(y,X) dy. (3.9) 

We denote the spherical harmonic coefficients of F, g’”en by the 

integral on the right hand side, by . Given F^ , we v;ish to 

solve for the coefficients of A . We have 

n 

<^1 ®0 + “0 ®2 = ''o ' '=1 = -“0 * 0^5 ' 


=2 ®1 + “S =3 = ' '=2 = -<2+f Co)A° , 


6° B° + a® bJ + C,., B® = F® , Cj 


-2 _0 
3 °0 ^2 ' 


B° + a°3 + C4 bJ = pO , C4 


2 nO 
" 5 °0 ^3 ' 


®k “k =2+2 + ^2 =2-2 = f 2 ' < i i ■’-2 - 


N-1 N-1 'N-1 ®N~3 'N-1 


» 


^ ®N “2 



The N+1 equations for F^, 0 £ k £ N+1, can be put into matrix 
vector form as a five-diagonal system in , 0 1 n £ N+1 : 



( 3 . 10 . 0 ) 

where 

°1“ "*^3 “0/^1 ^2 “ "'^4 ' 


The equations for P,'*' are 

K 

Cg b 1 + ai b 1 = pi , Cg = - , 

»2 ®2 + “2 = ^2 ' 

11 *11 11 2 1 
et b;: + at b. + c., Bt = Ft, c., = - c„ a- , 






We have a similar system for ; 

Cg b”^+ where Cg >» (-2 - Cg |) , 

®3^ * “3^ ®5^ + CgB-l=P;\ where . (-co 

^'k ®ki2 ®k-2 ” '■k^' " < >= < N-2 . 

»N-1 ®N-1 " "n-1 ®N-3 - ''n-1 ' ■ 

®N ■*■ ''^N ®'N-2 ■ ' 

All the other systems are of the form 



k=2, . , . , N , O.lO.k) 
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In O.lO.k) 


^ 4. nii 


K K * “k K■^2 - 4 


for 


I j I « 2,3f.,.f M-3 and k « j»j+l ; 


j I = N-2 and k = N-2 , 


4 4 + “k 4+2 -^k =k-2 - K 


ni 






for I j| *= 2,3,. . . , n- 4 and k - j+2 , . . . , N-6 ? 


4 4 ■= 4 


i jl = N-2 
for ^ 1 jl = N“1 
j|- N 


and k = N-1 » 

and k = N-1,N , 

and k = N ; 


4 4 + 4 4-2 ■= 4 


for 




Ijj = 2,-j, + . . , N-3 and k = N, N-1 , 


(jl = N-2 -and- k = N *. 


Each of the above linear systems are -penta diagonal. .They can be 
solved by the LU factorization algorithm given in Appendix A. 


After solving for the vi'e. wish to obtain the covariance 


function V * More specifically, we. are interested 

in the zonal covariances, , 


R(0,x) * r(0,XQ?e,XQ+ X) ; 


(3.11) 


I I i 'll lyiiiiit'iilf II 


r-r iimTiii t fr ii rr i A l ii i r i iii # n i d i i i i i ,' iik n' , 
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clearly R(0fX) is independent of X^. The zonal variances, 
R^(0,X), oottespondin, to t truncated aolution <■<«), are given by 


- 

Rjj (8,A) = I (e,y)4<”> (e,y +X)dy 


!_ .!. I j 

-n 


N n 
n»0 n*"0 m* 

'27T 
0 


;27T 

J git (m+m')y + m' X] 


dy 


N N 

^ I I n”^ n”ri ^“ilUX 


n=0 n'»0 m=-n 


I b; B- P™(u)P-(p) e- 


N 


« 2tt 

• 1 
n=0 

(B° 
' n 

P° 


V 





N 

N 

n 

bJJ b"”' 

n n 


2 1 
n=»l 

1 

n ' = 

•1 m==l 


„Li 


(3.12) 


Clearly, (3.12) does not provide a computationally efficient way 
of obtaining Rj^(0,X) for large N » 

In the following two sections we shall use the numerical approach 
of directly inverting (2.2) for fixed RHS, u = constant. R(6,X) 
will be computed by summing over correlation products of solutions 
for various « . A iast solver for the inversion of (2.2) with 
given RHS(w = const.) is presented next. 


X7. 

4 . A generalized fast Poisson solver 
A fast solver for the e<juation, 

+ c<e)]4. . f(e,A), (4.1) 

with C(0) « constant, is available in the NCAR software 
library (Swarztrauber and Sweet, 1975). This progrsOT was 
modified and extended to solve Eq.(4.1). The NCAR version takes the 
finite-difference approximation to Eq. (4.1) as 




- since. - i^e 


i ^ ^ _ 2(|). . + <J). + Ci}>4 . 

1^1 sinO^)^ ^0+1 ifj itO 1 iO 


f ( 0 / A j ) , 

A , j 


(4.2) 


where C - const , 

0^ = (i-l)A0 , A^ = (:)-l)AA , 

= «J>(0. , A.) , A0 = ir/M ,AA = 27 t/N , 

‘ l~l,2,...,M'^ 1 , J~1,2,...,N4‘1 . 
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In the NCM version, the option of solving (4.1) in a subdomain 
on the surface of the sphere exists. This option was removed in 
our version. Our extended program then takes C in eq. (4.2) to 
be Ci » C(0i). No futher modifications of the original program 
were required. The modified program was extensively tested 
(Table 1). 

We conclude from the numerical tests that the modified ver- 
sion of the solver has the desirable properties of the original 
one 5 second-order accuracy with respect to discretization error, 
and machine accuracy in linear system solving for a moderate-size 
mesh (32x64). 
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5. Monte“Carlo solution of the model equation 


We wish to solve Eg. (2,7) with the RUS P(e,iiw) h^ipg white noise, 

The solution is approximated by Monte-Carlo simulation. In other 
words, Eg. (4.1) is solved for an individual realization of the 
white noise process P(0,)k;w), namely f(0,A). A large number N 
of samples drawn from the process P(0,X;w) , 

Let <j)^(6,X) be the solution of (4.1) with f^(0,X) as rhs and let 

4 


(N) 

Our Monte-Carlo solution of (2.7) is ()> 

If {f^,n«l#2 , , . , ,n) are simply independent samples of P, 

(n5 

then a very large N is needed in order for ■ to be a good 

approximation to the true solution of (2.2); i.e., has 

a large variance around its expected value 4> , even for large N. 
In order to accelerate this convergence of to , we used 

a technique for variance reduction called stratified sampling 
(Appendix B) . 


Eg. (2.7) was discretized with 19 grid points in the meridional 
direction and 32 grid points in the longitudinal direction. The pole.^ 
are grid points, and the equator is a grid line. 

We obtained the zonal correlations from our simulation as 
follows. Let the normalized correlation r of an individual realization 
be defined as 

N N 

2it 


r^J0,A) 


4>^,(0,A')<|>^(e,X»+X)dX' 


tZir 2 

(}>,J0,X)dA 

K, 


(5.1a) 


0 
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notice the slight change of notation from (3*11). r^(OA) depends 
only On the separation in longitude, X, between two points (O^Xq) 
and (0 ,Xq +X) on a circle of latitude 0 , and not on their position 
Xq. This is due to the rotational invariance with respect to 
longitude of (2.2). Now the correlation function r(0,X) is approximately 

r'^^ (0,X) - sj I r . (Sab) 

Here we assume that the mean of <|> is zero. We took a large number N 
, of correlation realizations, r^^, and averaged them to obtain statisti- 
cally stable results. 

The Integrals in (5.1a) were evaluated by the trapezoidal rule, 


M-1 

Rj - R^.(e,jAX) » 4,^ AX, (5.2) 

Where we have dropped the argument 0 and the subscript nr, so that 
« <j) (0,kAX) . 

M-1 

The convolution product l be computed using 

Fast Fourier Transform ideas (Henrici, 1979) . Let 3R be the 
reversion operator? for a periK'‘?'ic sequence {x/} , it is defined 
as 

(3Rx)^ = . (5.3a) 

Let 3F be the discrete Fourier operator defined, for two sequences 
{x^> and of period M by 3F x = & ? component-wise 


0 

m 


M 


M-1 , 

I X 


k 


t 




with w *• exp(^~ ) . It follows that the inverse Fourier 
is defined by l » MK R - M E K . 

Consider the Hadamard product x*y of two sequences x 
X • y * {xjj yjj}f and let 


U“T x,V«r y 


■* r*- 


The Fourier transform of x*y is then 
(* <*•?>>» - H *K “ 


-km 


, M“1 

M-1 

M I X. 

” k*0 ^ 

1 ^SL 

M~1 

Z„ Vj 
Jl»0 

> M®*! 
i. f - 

“ k-o 

M-1 


1 

f*0 

Vn . 


£k . “km 


-(m-i)k 


For any two M-periodic sequences U « {Uj^} and V » » 

the cont'olution product C 


M“1 M-1 

=- = L \ “m-k = J. ''k ''m-k 


m 


k»0 


k=0 


We denote this convolution product by 5 = 0 * ^ . 
Thus 

**V' 

t (x*y) » I X * Ey . 
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operator 
and y , 


we define 


( 5 . 4 ) 


n* 


Taking into account that ^CQ#X) ia 2»f -periodic with 
reapect to K, we have that 


M"3, M~1 M- 

1,^0 *>= “ do *>= ■ do 


or 

»rV 

R » # * EiJ) 

- 3c ( $) ) 

«* JE((E E?) ’ (E E E J) ) 

« E E ( (E E$ ) • (eJ) ) 

« M E“*( (E e|) ’ (e|)) 

«ME"^{(r|)** (e|)) , 


After calculating the correlations rj. r r (0/ jAX) for each 

J K 

latitude D at uniform angular distance intervals hX , we use spline 
interpolation to obtain the correlations 'r. (e,s) in terms of soherical 

1C 

distance s around a circle of latitude, at regular intervals Ac® 200km. 

Our experiments showed that Monte-Carlo simulation/ even 
using variance reduction/ converged very slowly. It required a 
large number of realizations. But the solution had the expected 
behavior, namely R(0,s/a) - exp(“?/SQ), where SQC0)~sinO , 
and a is the radius of the earth. 
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® ’ An equation for the covarianc® function 

Given the linear equation with deterministic coefficients 
(2.7) for ^ it is easy to obtain a linear equation for 

tho covariance function ' defined in (2.3) . Writing 

(2.7) symbolically as 


L(Ji » F(5 jw) 


(6 . la) 


with 


L - (y^-c^sin^e) , 


(6.1b) 


we are Interested in the ensemble average of ^ 4> (§ 2 • 

Let Li, L 2 be the operator L written with respect of the 
position vectors and ? 2 # 

and Pj^^2 * ^^^1,2*“^ ' 

*^2 ^1 *^1 “ *^2 ^1 ' 


applying L 2 taking the ensemble average, we have 

E L 2 4»2 % ‘*’1 “ ^2 ^1 ^2 


E L 2 4)2 = E F 2 s 6(5i“ § 2 ^ 


Here we used the fact that E, L^, L 2 are all linear operators 
and commute with each other, since they operate on the independent 
variables w, and ^ 2 ' respectively. We thus obtain for 

r(5^, ?2^ == ®^i ‘^2 




W,£33.ag-.fetra^. 
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* 

the deterministic linear equation ! 

A rigorous derivation of (6,2), under suitable assumptions on the 
operator L and data P, appears in Bdcus and Cozzarelli (1976). 

We are only interested in computing 

R(0,X) = r(0,O;0,X) . (6.3) 

To solve (6.2) for R(0,X), we start by solving numerically 

for each First, the solution H(?^; for given 

obtained for Xg ~ More precisely, for = (k2A0 ,O), 

k 2 ” 0,1,..., K“l, we compute H(kj_A9, ^ 2 ) ' *^1 0,1,..., K-1, 

~ 0,1,..., J“1 • 

i 

The solution other ^2 “ j2^^) ^ 

j2 = 1,2,. . . ,J-1, k 2 = 0 , 1 ,..., K-1 is then obtained from the ^ 

previously computed H(5x'?2^ 

H(kxA0, j^AX; k2A0, j 2 AX) = H (k^AO , ( j x" j 2 ^ ? k2A0,O) . (6.5) 

In other words, only K "inversions" of Lx are needed, rather 
than KJ in solving (6.4), viz., O(K^J) operations rather than 
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0(k2j 2) have to be performed. Furthermore, (6.3) shows that 
we need only solve 

L2r(5i; 52) - h(?i:52) <6.6) 

for H(8,0;52)» since we ar<i only interested in R(9,x), rather 
than in the full F(l5i;€2). Hence, only K "inversions" of L2 
are needed, leading to a total 0 (k2j) operation count to obtain 
R(0,X) from (6.2). In our numerical solution of (6.2) directly 
is much more efficient in order to compute R(6,X) than a Monte- 
Carlo solution of (2.7). First, N, the number of realizations 
necessary for a good approximation (j)<N) to is considerably 
larger than 2K. Second, here the covariances are obtained directly, 
without the need for computing convolution products of (j)(N) ; 
the later requires 0(Kj2iog J) additional operations. 

The numerical results we present in the next section depend 
only on the form of Eq. (2.7), not on the method used for obtaining 
it solution and the covariance of this solution. The results were 
obtained by the most efficient method, that described in the pre- 
sent section. Some of the results were futher confirmed by using 
Monte-Carlo simulation with variance reduction (Section 5 and Appen- 
dix B) . The zonal covariances R(0,X) obtained by solving (6.2) 
were divided by the local variance to obtain zonal correlations 
r( 0, X) as in (5.1) . 
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7. N utnerlcal results 

i* 

The starting point of this study was the development and 
testing of a model which will reproduce the latitude dependence 
of s5onal correlations of observed in GHA. The model is 

governed by Eqs. (2.7) and (6.2). 

^ • Data sets 

The model was tested first using the data on which Figs. 2a-d 
of GHA were based, namely vertically averaged correlations of 
the difference field, TO - T^ , between satellite-derived and 
model forecast temperatures. This test was based on the hydro- 
static connection between temperature and layer thickness, 
which should make and vertically averaged T interchangeable 
at the level of approximation of our model. Model correlations 
exhibited the same tightening with increasing latitude as that 
of the data, for any reasonable value of the parameter Cq. 

After this preliminary test, we proceeded to study If our 
model results would fit observed 500mb geopotential height 
data. We considered first NMC analyses of 500mb heights for 
the DST-6 period, January-March 1076, and for January 1978, 
available on tape at GLAS. Zonal correlations for the analyzed 
fields every 12h over these periods were computed at different 
latitudes and averaged over time. The averaged correlation 
curves showed the typical damped cosine (Gandin, 1963; Thiobaux, 
1977) or Bessel function behavior (Jq: Rutherford, 1972) of 
atmospheric fields. In particular, the curves crossed the 
s-axis, becoming negative at some distance sj. =sx ( 9 ) 
which decreased with latitude. It is clear that ou^ ^ h 
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not simulate such behavior: the covariance function 
given by (6.2) cannot become negative. Indeed, the operator 
(-L) in Eq. (6.1) is monotone , i.e., -L (|> ^ 0 implies <j> 
this property is closely related to L satisfying a maximum 
principle . Thus ^(5i-?2) Z ^ implies -L2r(?3^,52) >. 0 
and hence r( 55 ^, 53 ) 2 i.e., r is nonnegative for all 

values of 5i and 52* More rigorous proofs, not involving S - 
functions, can be given, using the fact that r(|i,?2) is 
closely related to the Green's function of the monotone operator 
Ll L2. 

The physical reason that correlations of observed geopoten- 
tial data become negative at distances larger than si is the 
presence of the planetary vorticity in the potential vorticity 
Eq. (2.1). The planetary vorticity varies with latitude and 
it is the v9f/80 term that gives rise to Rossby waves, which 
have a characteristic horizontal wavelength of the order of sev- 
eral thousand kilometers. These waves dominate large-scale, 
mid-latitude flow and hence the statistical structure of mass 
fields at these latitudes. The variable planetary vorticity, 
however, is not present in the geopotential error equations 
(2.3) or (2.5). This makes it physically plausible that the 
horizontal error correlation should never become negative in 
our model. 

Based on the evidence of the correlations for the difference 
fields T° - T^ , we turned our attention to the model forecast 
fields of 500mb heights. One of the assimilation experiments 
performed at GLAS with DST-6 data, namely the one which, in 
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addition to conventional data, assimilated tempera tut'e data 
from two satellites by a time-continuous statistical method 
(experiment S2a, cC. GHA, Table 1) was used, Thtee-day forecasts 
had been performed from initial states produced by the assimila- 
tion cycle every 48h from OOOOGMT 1 February 1976 till OOOOGMT 
4 March 1976. Table 2 shows the availability of forecast 
fields at synoptic epochs for all the forecasts, at the time 
our computations were carried out. 

b. Model validation 

We computed correlations for the difference fields ij)° - 
using NMC analyses as and our experimental forecast fields 
as *1) for 12h 24h, 36h, and 48h forecasts (Table 3). These 
correlations were then compared with our model correlations, 
using different values of the nondimensional parameter Cq. 

As in the case of the temperature difference fields, model 
results and data were in good agreement for a certain range of 
values of Cq* The result which is apparently best for all 
error measures is boxed. This value seems to increase slightly 
with forecasting interval: it is ilO for 12h, 150 for 24h, 170 
for 36h and 190 for 48h. For the combined 24h and 36h data 
set, it lies between 160 and 170. 

c . Effects of vertical stratification 

The stochastic model (2.7) gives good qualitative agreement 
with zonal correlations of experi nental forecast errors, provided 
we use values of the constant Cq oetween 100 and 200. From 
Section 2, 

= 4 R 2^2/gp^ (7.1) 








29 


For the external mode of the atmosphere, D 10 km, this expres- 
sion takes a value of cq - 10. The lack of agreement with our 
results indicates that external, barotropic motions do not 
dominate the error field growth. 

As a first step in accounting for the effect of vertical 
variations, we can use in (7.1) a reduced value of g wity, 
g'«Sg, that takes into account the stable stratification S, 

S=(D/0s) dQs/dz, 

of a standard atmosphere with potential temperature Qg = Qs(z) 
(Pedlosky, 1979, Sec. 6.5). A typical value of the stratifica- 
tion parameter is S*=0.1, which corresponds to a Brunt-Vaisala 
frequency of N s 10“^ sec"^, where N^-(g/0g)d0g/dzi. 

Redefining Cq as: 

Cq = 4n2a2/g'D, (7.2) 

and substituting a value of g' based on S=0.1, we obtain c^ s lOO. 
This is in much better agreement with our results. 

A rigorous analysis of stratification effects would require 
a full three-dimensional separation of variables. Such a devel- 
opment does not seem to be warranted by the poor accuracy of the 
data. Stratification of the data by height would also impose 
serious restrictions on their statistical reliability, due to 
limited sample size. 

We can interpret our results as indicating that forecast 
error growth is dominated by baroclinic motions in the atmosphere. 
Baroclinic instability is the most importa^nt dynamical instabi- 
lity in the extratropical atmosphere. Hence this aspect of 
our results is not entirely unexpected. 


30 . 


8. Discussion and conclusions 

We have derived a stochastic-dynamic model for the global 
structure of the atmospheric mass field forecast error. The 

model is governed by a stochastically perturbed Helmholtz-type 
equation (2.7). The covariance function of the model's solutions 
has been shown to be governed by Eq. (6.2). The study originated 
in the observation (GHA) that zonal correlations of mass fields 
exhibit a strong latitude dependence, with radii of equal corre- 
lation becoming smaller at higher latitudes. An analysis of 
the potential vorticity equation, both for the atmosphere and 
for numerical weather prediction models, led to the derivation 
of our stochastic-dynamic model for the error field. The model 
supports the heuristic interpretation of the tightening of the 
correlations with latitude given in GHA, namely that the typical 
correlation radii vary with latitude in the same sense as the 
Rossby radius of deformation. 

The observation (GHA) that correlation fields did not ex- 
hibit a strong vertical dependence led us to choose a single 
model parameter to represent the vertical structure. The value 
of this parameter, Cq, was determined so as to give the best 
quantitative fit between the second-order statistics of the 
model and of the data. The empirically determined value is 
consistent with baroclinic motion dominating the error growth 
field. The barotropic mode does not seem to play a major role. 

The model results are in good agreement with actual numeri- 
cal weather prediction errors, both for temperature and geopo- 
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tential fields. The differences between model and data are 
typically smaller than those between the data of one hemisphere 
and those of the other. We intend therefore to use the global 
correlation functions produced by the present model in the 
development of the GLAS statistical assimilation method. This 
should eliminate the inconsistency produced at present by the 
use of ad-hoc , meridionally stratified, empirical correlation 
functions (GHA). 

The agreement between our stochastic model and actual fore- 
cast errors fields for 12 to 36 hour periods validates the as- 
sumptions on which the model was derived. Within this period, 
the difference between the potential vorticity fields of the 
atmosphere and of the numerical model forecast used in the 
comparison is well represented by white noise. 

For periods shorter than 12 hours, the lack of balance in 
the initial data generates fast inertia-gravity waves in the 
model, which violate the quasi-geostrophic assumption implied 
in the potential vorticity equation. The use of a strongly 
damping time scheme, the Euler-backward scheme, in the GIAS 
second-order mode3, makes these waves negligible after 12 hours. 
The use of an effective initialization scheme could eliminate 
this restriction. 

For periods longer than 36 hours, the forecast errors become 
so large that the model equation (2.2) ceases to be a good re- 
presentation of the atmospheric governing equation (2.1). The 
length of this limiting period is obviously model dependent. 

The forecast model used in obtaining our forecast error data, 
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the GLAS 9*^level 4® latitude by 5" longitude, second order 
model, has rather large truncation errors. With higher resolu- 
tion or higher accuracy models (e.g, Kalnay-Rivas and Hoitsma, 
1979) we may expect that the stochastic model will remain valid 
for longer periods. 

The observation that our stochastic model based on white- 
noise forcing fits Southern Hemisphere errors better than those 
of the Northern Hemisphere is also an interesting result. It 
indicates the validity of an increasingly well accepted points 
that the lack of proper parameterization of the planetary- 
scale forcing can have a very important effect on the forecast 
errors . 

The interplay between atmospheric statistics and dynamics 
which is stressed by this work points the way to further studies 
along similar lines. Allowing for such interplay in modeling 
efforts might help to improve our knowledge of both deterministic 
and stochastic aspects of atmospheric behavior. 

A concrete step along this road would be to investigate 
the stochastically perturbed potential vorticity equation (1.1) 
itself, rather than its steady-state form (2.2). This would 
be a nonlinear, time-dependent Langevin-type equation for 
large-scale atmospheric flows. Its study appears considerably 
more difficult than that of the present model, but still acces- 
sible by Monte-Carlo simulation with variance reduction. We 
hope to report on such a study in a future publication. 
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Appendix A . Algorithm for penta“diatfonal system 


We solve 

Ax * b (A.l) 

where A is penta-diagonal, by LU factorization 

A » LU . 

The factorization can be derived from the well-known one for 
tri-diagonal systems, by noticing that actually cdd and even 
variables in (A.l) decouple. Hence . 



The algorithm can be performed in the following three steps. 
S tep 1 » ^ 1 ^ ^2 ” ^ 2 ^ 

^\-?/‘^k-2^ ®k~2 ^ 
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Stop 2 (Forward elimination# i.e.# solve Ly * b) ; 

yi - bj^ , yj - b , ■ 

” ‘’k " \-2 ^k-2^S-2 ' • 

Step 3 (Back substitution# i.e»# solve Ux ■» y) i 

><h - l'n''=n ' Xn-1 ” yn-l/=n-l ' 

** ~ “k ^k+2^^^k ' ^ “ n-2,n~3# , , • #1 . 


I 


r 


L 
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Appendix B . Stratified samplinq 

Monte-Carlo simulation is a method for the solution of 
stochastic model equations using sample realizations. The 
individual realizations are computed by using pseudo-random 
number generators to specify the random functions prescribed 
as data of the equation. Basic random number generators produce 
variables which are to a good approximation uniformly distributed 
on the interval [-1,11 and stochastically independent. 

We used the following stratified sampling procedure 
(Kleijnen, 1974, Ch.3), bet N be the number of realizations in our 

t 

sample. We classify the random value taken by any sampled RHS, 
fj (0/X), at each grid point (0^,1,.) into exactly one of N classco. 

A ■•JR 

These classes are formed by dividing the range [-1,11 into N non- 
overlapping exhaustive intervals, we picked the range of each class 
to be of equal length. The functions fj(0,X) and fj^(0,X) are made 
dependent for j k, j, ks{l,2,,.., N) , by taking at each grid 
point fj(0,X) and fj^(0,X) from different classes. The objective 
of stratifi.cation is to lead to variance reduction. 

We solve Eg. (4.1) with f(0,X) « f^(0,X) to obtain 
(}>., i « 1,2,..., N. Thus the p. themselves will be dependent. 
After experiments with both independent and stratified sampling, 
the superiority of the latter was clearly established. A sample 
size of N 200 gave satisfactory results when using stratified 
sampling; independent sampling with N ~ 700 gave results whose 
variance was still much too large to be satisfactory. 
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Table 1. Test of: our extended fast Poisson solver. 


Trial 
• function 
t|)(0,X) 

Form of 
forcing 
function 

C 

1 

M 

1 

N 

1 

Max 1 (() 1 

Max| 4 ,| 

» 

2 

cos 0 cosX 

analytic 

- sin^0 

16 

32 

1 

8.16 

-2 

II 

II 

II 

32 

64 

1 

2.3 

-2 

II 

II 

-3.6 +7 

16 

32 

1 

8.3 

-8 

II 

II 

-10.0 

16 

32 

1 

3.57 

-2 

II 

finite 

difference 

-sin^e 

32 

64 

1 

2.9 

-12 

2 

-cos 0 cos 3X 
48 

analytic 

II 

64 

128 

0.02 

1 

7.19 

-3 

cos4> .cos 0 

II 

2 

- 10 sin 0 

32 

64 

1 

1 

1.18 

-2 

II 

II 

-5 sin^e 

32 

64 

1 i 

1.36 

-2 

II 

II 

• 2. 
-sin 0 

32 

64 

1 

2.3 

-2 

(sin 0 sinX • 
2 

cos 0 cosX) 

II 

" 

32 

64 

< 0.25 

9.14 

-4 


The RHS of (4.2) corresponding to a given trial function was 
computed either as grid point values of the residual of that 
function (4.1) ("analytically") or as the residual of that func- 
tion in (4.2) ("finite-difference"). M is the number of grid 
intervals into which the meridional coordinate 6 was divided; 

N is the number of grid points in the longitudinal direction. 

A number a.b + c means (a.b) x 10 
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Table 2. Forecast fields from the DST-6 experiment S2a (GHA) 
which were available for our computations. 


Date 12h 24h 36h 48h cOh 72h 



Forecasts were initiated at 0000 GMT on the day indicated in the 
first column. They were carried out for 72h each, but some of 
the existing prognostic fields were lost in data storage, trans- 
mission or retrieval. The columns in the table indicate the 
availability of the fields at successive synoptic times for each 
forecast. NMC objective analysis were available at all of the 
corresponding synoptic epochs. Beyong 48h there were too few 
fields available to constitute a statistically valid sample. 




Table 3. Comparison between model correlations and 
correlations of the data fields. 
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c 

o 

(Ar) 


Li (Ar) 

(-2 ^Ar) 



12.15 

10 

.322 

.379 

.379 

.442 

.84 9 

.857 

12.14 

50 

.132 

.176 

.207 

.257 

.640 

.685 

12.13 

70 

.108 

.143 

.183 

.221 

.589 

.647 

12.12 

80 

.101 

.134 

.178 

.212 

.569 

.633 

12.11 

90 

.0967 

.129 

.175* 

.207 

.552 

.620 

12.10 

100 

, .0943 

.126 

.176 

.206* 

.536 

.609 

12.9 

1 110 

1 .0938* 

.125* 

.180 

.208 

.522 

.599 

12.8 

120 

.0946 

.125* 

.185 

.213 

.510 

.590 

12.7 

14 0 

.0981 

.130 

.201 

.226 

• .488 

.575 

12.6 

160 

.103 

.136 

.219 

.24 2 

.470 

.562 

12.5 

170 

.106 

.140 

.229 

.252 

.462 

.556 

12.4 

180 

. 108 

.144 

.239 

.261 

.455 

.550 

12.1 

190 

.111 

.147 

.248 

.270 

.448 

.545 

12.2 

230 

.122 

.162 

.289 

.307 

.424 

.528 

12.3 

270 

.133 

.176 

.329 

.342 

.405 

.515 

24,8 

140 

.074 5 

. 0986 

.153* 

.172 

.488 

.575 

24.7 

150 

.0744* 

. 0970 

.155 

.171* 

. 479 


24.6 

160 

. 0748 

.0965* 

.159 

.172 

.470 

• iij U O 

.562 

24.5 

170 

.0756 

.0968 

.163 

.174 

.462 

.556 

24.4 

180 

. 0768 

.0978 

.169 

.178 

.455 

.550 

24.1 

190 

. 0784 

.0994 

.174 

.182 

.448 

.545 

24.2 

230 

. 0860 

.109 

.203 

.206 

.424 

.528 

24, 3 

270 

. 0939 

.120 

.232 

.234 

.405 

.515 



Table 3. (Continued) 


43 




(Ar) 




) 

36.6 

160 

.0898 

% 

.114 

.191* 

.203 

• .470 

.562 

36.5 


.0891* 

.112 

» 

.193 

.202* 

.462 

.556 

36.4 

180 

.0888 

.112* 

.195 

.203 

.455 

.550 

36.1 

190 

.0889 

. 112 

.199 

.205 

.448 

.545 

36.2 

230 

.0918 

• X. 5 

.217 

.218 

.424 

.528 

36.3 

270 

.0965 

.123 

.239 

.238 

.405 

.515 

48.6 

160 

.0952 

.123 

.202* 

.219 

.470 

.562 

48.5 

17 0 

.0937 

.120 

.203 

.216 

.462 

.556 

48.4 

180 

,0928 

.117 

.204 

.213 

.455 

.550 

48.1 

FTFol 

.0923* 

.116 

.206 

.212* 

.448 

.545 

4 8.2 

230 

.0928 

.113* 

.219 

.214 

.424 

.52 8 

A Q 0 

270 

.0951 

*11^ 

.235 

.226 

.405 

.515 

24/36.8 

140 

. 0808 

.104 

.166 

.182 

.488 

.575 

24/36.7 

150 

.0792 

r 

.102 

.165* 

.179 

.479 

.568 

24/36.6 

160 

.0786 

.0999 

.167 

.178* 

.470 

.562 

24/36.5 

[Tfo^ 

.0785* 

.0991* 

.17 0 

.178* 

.462 

.556 

24/36.4 

180 

.0786 

.0991* 

.173 

.180 

.455 

.550 

24/36.1 

190 

.0791 

.0996 

.177 

.183 

.448 

.545 

24/36.2 

230 

.0830 

,106 

.196 

.200 

.424 

.528 

24 36.3 

270 

.0887 

.115 

.219 

.224 ■ 

.405 

.515 

The first column 

indicates 

the label 

of the numerical 

experiment, 

based 


on the data set and on the value of Cq used. The data sets are identified 
by the length of the forecasts on which they were based; 12h, 24h, 36h, 
48h or combined 24h and 36h. The second column gives the value of Cq. 

The next four columns contain the various norms of the difference Ar 
between model and data correlations (Eqs. (7.1) and (7.2)), These columns 
contain a star next to the minimum value of each error measure fo^ 
corresponding column and data set. The values of Cr, which co>'’*os-'onc: to 
the largest number of stars in their own row and in the two adjr-cent rows 
are boxed. The last two columns give the norm.s of the model correlations 
themselvp';. 
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Fig. Id 















